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ABSTRACT: 

Molecular dynamics simulation was used to estimate the interfacial thermal 
(Kapitza) resistance between nanoparticles and amorphous and crystalline polymer 
matrices. Bulk thermal conductivities of the nanocomposites were then estimated using 
an established effective medium approach. To study functionalization, oligomeric 
ethylene- vinyl alcohol copolymers were chemically bonded to a single wall carbon 
nanotube. The results, in a poly(ethylene- vinyl acetate) matrix, are similar to those 
obtained previously for grafted linear hydrocarbon chains. To study the effect of non- 
covalent functionalization, two types of polyethylene matrices. — aligned (extended-chain 
crystalline) vs. amorphous (random coils) were modeled. Both matrices produced the 
same interfacial thermal resistance values. Finally, functionalization of edges and faces 
of platelike graphite nanoparticles was found to be only modestly effective in reducing 
the interfacial thermal resistance and improving the composite thermal conductivity. 


1. Introduction 

Polymers with improved thermal and electrical conductivity would have many 
applications [1]. Although addition of highly-conducting nanoparticles to polymer 
matrices is an active area of research, it is often found that thermal conductivities of the 
resulting nanocomposites fall well below a rule-of-mixtures-type additivity [2], The 
reason for this is assumed to be a large interfacial (Kapitza) thermal resistance between 
the nanoparticle and the surrounding polymer matrix. Modeling has been used to explore 
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trends in this resistance with particle shape, particle size [3-5] and interfacial molecular 
structure [6], 

Carbon nanotubes (CNT) are known to have extremely high thermal conductivity 
[7], and research in measurement and modeling of carbon nanotube composite materials 
is currently very active. In a previous paper [8], functionalization of single- walled 
carbon nanotubes (SWCNT) by covalent grafting of alkane chains to the tube walls was 
studied. It was shown using molecular dynamics (MD) simulations that functionalization 
could reduce the interfacial thermal resistance with a vinyl polymer by an order of 
magnitude. 

The present paper applies similar methodology to a new grafting architecture. In 
addition, because of the difficulty of controlling covalent functionalization and its 
potential for damaging the inherent properties of the CNT, possible non-covalent 
interactions are also explored - specifically, the effect of ordering of chains at the 
interface (coiled vs. extended-chain crystalline). Finally, a less-expensive filler, graphite 
nanoplatelets (GNP) is considered [9, 10], 

The multiscale approach involves atomistic MD simulations to establish key 
structure-property relationships for interfacial thermal resistance. To do this, the grafting 
density, < 7 , and length, n, of linear hydrocarbon chains (-C„H 2 „+i) covalently bonded to 
the nanoparticle are systematically varied and, in a separate study, the matrix ordering is 
specified. Effects on the thermal conductivity of the composite are then predicted 
parametrically using an effective medium approach. Key parameters in these analytical 
predictions are nanoparticle aspect ratio, volume fraction, and interfacial thermal 
resistance. 
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In Section 2, the methods for constructing the models and performing the analyses 
are described. In Section 3, the results from these analyses are presented and discussed. 
Section 4 summarizes the conclusions. 

2. Method 

The current work involves three steps. First, atomistic structures representative of 
the functionalized nanoparticle in the polymer matrix are generated using standard MD 
techniques. Next, the methodology in references [3-5] is applied to calculate the 
interfacial thermal resistance, Rk, in these structures. Finally, these Rk values are used in 
an analytical model of the corresponding nanocomposites. 

2.1. Structure preparation 

Three sets of structures were prepared. Due to their complex architectures, care 
must be taken in order to produce results that are representative of the actual structure. 
The first set involves ethylene vinyl alcohol (EVOH) chains grafted to a single-walled 
carbon nano tube (SWCNT) embedded in an ethylene vinyl acetate (EVA) matrix. Each 
simulation cell was constructed as follows: 

A (6,6) SWCNT consisting of 577 carbon atoms was functionalized with 0, 5, or 
10 carboxylic acid groups (corresponding to grafting densities < 7 = 0, 0.00338, or 
0.006777 A" 2 ). Ten EVOH copolymer chains, each consisting of a random sequence of 
10 vinyl alcohol monomers and 15 ethylene monomers were then generated. To model 
the grafting chemistry, a short constant volume, constant temperature (NVT) MD 
simulation was run in order to allow the copolymer chains to condense around the 
SWCNT. Then covalent ester links were formed between vinyl alcohol groups on the 
copolymers and the carboxylic acid groups on the SWCNT as indicated by Eq. (1). 
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Hydroxyl groups are in excess in this reaction, and all carboxylic groups react; water 
molecules were removed from the simulation. This functionalization architecture is 
different from the end-grafting of linear chains that was modeled in ref. [8], 

SWNCT-COOH + HO-Copolymer -> SWCNT-CO-O-Copolymer + H 2 0 (1) 

The functionalized tube was embedded in an amorphous ethylene- vinyl acetate 
(EVA) copolymer matrix consisting of 60 chains of EVA containing 25 repeat units each. 
The 15% vinyl acetate monomers by weight were inserted randomly. Since the vinyl 
acetate mole fraction was only 5.5%, each chain contained only 1-2 vinyl acetate 
monomers. Each of the three grafting densities was simulated 3 times, for a total of 9 of 
these structures. Bulk density was effected through application of periodic boundary 
conditions [11]. In all the simulations, the PCFF force field [12] was used with the 
LAMMPS [13, 14] MD simulation software. The SWCNT was considered to be 
continuous through the z-axis periodic boundary coordinate; its atoms were not otherwise 
constrained in any way. 

The bulk atomistic nanocomposite model was constructed by compression from 
low density. This procedure began by placing the SWCNT with its grafted chains and the 
polymer molecules in a box with periodic boundary conditions. Initially, the density was 
about 5% of the bulk density of 0.8 g/cm 3 . To prevent the polymer chains from 
collapsing, the nonbond pair potential used to describe the forces between atoms which 
are not chemically bonded to each other was scaled by a factor of 0.001 and a short cutoff 
(3.0 A) was applied with electrostatic interactions between the partial charges on each 
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atom turned off. To allow the polymer chains to diffuse through the periodic box and 
achieve random configurations, the MD simulations were run for 50 ps using a 1 fs 
timestep at 500 K. The box volume was reduced linearly with time to give a density 
about 10% of the final bulk density. A second MD run condensed the system to a density 
of about 50% of the final bulk density over 50 ps at 500 K. During this run, the unsealed 
nonbond pair potentials were then applied with a cutoff of 9.0 A, and electrostatic 
interactions between the partial charges on each atom were included via the Ewald 
summation technique. At the end of this run, energy minimization was performed on the 
system, followed by constant pressure MD simulation at 300 K and 1 atm of pressure. 
After that, MD was applied for 10 ps while increasing the pressure to 100 atm. This was 
followed by 50 ps of MD simulation. The pressure was reduced from 100 to 1 atm over 
10 ps, followed by a MD run at 1 atm and 300 K for 100 ps. Such equilibration 
procedures are typical for the simulation of polymers[15, 16]. Running MD at high 
temperatures usually allows the polymers to relax out of any nonphysical configurations. 
Running MD with elevated pressures helps to achieve an appropriate condensed phase 
density. 

The second set of structures involves a SWCNT surrounded by two different 
polyethylene (PE) matrices. The first of these, amorphous polyethylene, was constructed 
in a manner identical to the one described above but using 78 chains of 45 carbon atoms 
(C45H92) as the matrix. A second model was constructed in which the polyethylene 
chains were aligned parallel with the SWCNT axis. The chain length (45 carbon atoms, 
C45H92) was chosen to fit in the space available, with the SWCNT extended through the 
z-axis periodic boundary conditions. The 78 C45H92 chains were arranged parallel with an 
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initial density of about 70% of the final equilibrium density. Energy minimization with a 
convergence criterion of 10 kcal/mol was first applied to the configuration. MD was 
applied for 10 ps while increasing the pressure from 1 to 100 atm at 300 K. This was 
followed by a compression at 100 atm for 50 ps at 300 K. The pressure was reduced 
from 100 to 1 atm over 10 ps, followed by an MD run at 1 atm for 400 ps and 300 K. Fig. 
1 shows axial views of these two polyethylene-SWCNT atomistic models. The ordered 
matrix shows the polymer chains retaining their orientation parallel to the axis of the 
nanotube after the equilibration run. 

The third set of structures involves the GNP in an EVA matrix. Because the 
planar dimensions of a GNP are still very large compared to the typical volume of an 
atomistic MD simulation, the platelets are simulated as (infinite) sheets extended through 
the periodic boundary conditions. Two surfaces (basal plane and edge) are simulated 
separately. To model the flat surface, 3 sheets of graphite each containing 836 carbon 
atoms were oriented in the xz plane, with periodic boundary conditions along the x and z 
directions. The normal to the graphite plane is the y coordinate as shown in Fig. 2. To 
model the edge of the graphite particle, 13 sheets of graphite, each containing 228 carbon 
atoms, were oriented in the xz plane with periodic boundary conditions along the x and y 
directions as shown in Fig. 3. 

Surfaces and edges are functionalized by forming covalent chemical bonds with 
short linear hydrocarbon chains (GNP-C„H 2 „+i) with a range of grafting densities. The 
GNP/EVA nanocomposite models are prepared via the same methodology used for the 
SWCNT nanocomposites, with the exception that the graphite atoms are held fixed 
during the condensation procedure. The final model uses periodic boundary conditions in 
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the xyz directions. The polymer is sandwiched between the graphite structure in a 
condensed phase. 

2.2. Simulation of the heat transfer 

Following the last MD run of the atomistic structures at 300 K described in 
Section 2.1, Rk, the interfacial thermal resistance was calculated [3-5] by instantaneously 
heating the nanoparticle atoms to a temperature of 500 K and then monitoring the 
difference in temperature, AT, between the nanoparticle atoms and the matrix atoms 
during a constant energy simulation (NVE ensemble). A plot of ln(A7) vs. time yields a 
slope which is the negative inverse of the characteristic decay time, T. Eq. (2) then gives 
the interfacial thermal resistance, R K , where c T /A T is the heat capacity per area of 
nanoparticle surface. This ratio (cr/Af) is calculated assuming a value of 0.71 J/gK for 
cr[17]. This results in a value of 5.6xl0~ 4 J/m 2 K for cj/Aj for both the SWCNT and the 
flat surface of the GNP while a value of 3.1xl0~ 4 is obtained for the edge surface of the 
GNP. 


Rk 


T 

( c T / Aj. ) 


( 2 ) 


2.3. Analytical modeling of the thermal conductivity 

As was done previously for functionalized SWCNT nanocomposites [8], thermal 
conductivities of nanocomposites were calculated using the interfacial resistance values 
from the heat transfer simulations. For a dilute dispersion of randomly oriented prolate 
ellipsoidal inclusions (dimensions a\= a^ > af), an effective medium[18, 19] approach 
yields Eq. (3) 

K =K 3+/[2fl,(i-A l >+/g»a-.i„>] 

' " 3-/PAA, + A,i»] 
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where K e is the thermal conductivity of the composite, K m is the thermal conductivity of 
the matrix and/is the volume fraction of the inclusion. The indices T’ and ‘3’ refer to 
the dimensions a\ and a?, respectively. Fig. 4 depicts this ellipsoid. The subscripts (ii) in 
La refer to the coordinates of the GNP, where ‘3’ is the direction perpendicular to the face 
of the ellipsoid shaped GNP particle and fin is defined as: 


P« = 


EiZlrn 

K m +L U (K‘-K m ) 


(4) 


and La is defined as: 


T) T) 

L, = L,, = , + , w cos 1 p, 

Ml 2(p 2 -1) 2(l-p 2 ) 3/2 

(5) 

^33 ~ 1 _ 7L, 1 

(6) 


respectively, where p=apa\ is the aspect ratio of the ellipsoid. The values K c u , are 
defined as 


K. 


K n=- 


2R 


1 + 


parallel 

~K C 

parallel parallel 

a , 


(7) 


K. 




c perpendicular 


2R 


K 




perpendicular c perpendicular 


( 8 ) 


where K r and K r are the thermal conductivities of the graphite layer in the 

perpendicular and parallel directions respectively. In evaluating Eqs. (2)-(8), we used 
K m = 0.2 W/mK, K = 390 W/mK [17], and K =2 W/mK which are the 

c parallel c perpendicular 

values given for pyrolytic graphite. This should be a reasonable estimation for the GNP 
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are those 


properties. The interfacial thermal resistance values R„ and R K 

A perpendicular A parallel 

calculated from atomistic MD simulations. 

3. Results and Discussion 

The following three subsections describe results for the EVOH-functionalized 
nanotube, the functionalized nanotube in ordered and amorphous polyethylene matrices, 
and GNP nanocomposite modeling, respectively. 

3.1. EVOH functionalized nanotube/EVA 

The nanocomposite model was described in Section 2.1. Fig. 5 shows the 
interfacial thermal resistance, Rk, as a function of grafting density, a. Error bars shown 
are standard deviations among runs for all 3 independent configurations. Although the 
effect of ester- linked EVA on Rk is substantial, it is noticeably less than the effect seen 
with directly grafted linear hydrocarbon chains [8], Thus as might have been expected, 
the coupling of SWCNT vibrations to the polymeric matrix [5] is seen to be sensitive to 
the chemical bond topology adjacent to the tube. These results also serve to illustrate the 
value of MD simulation in providing guidance on chemical functionalization with the 
greatest potential payoff. 

3.2. SWCNT/ Amorphous and Ordered Polyethylene (PE) 

Two different types of PE/CNT composites were constructed. In one type, the 
matrix chains were amorphous random coils; in the other, the PE chains were oriented 
parallel with the axis of the carbon nanotube. Fig. 1 shows views along the z-axis of the 
amorphous and ordered matrix atom configurations. As before, there were 5 thermal MD 
runs for each of 3 independent configurations. The values obtained for the two matrix 
types are identical (10.9 ± 0.63 xlO" 8 m 2 K/W for the amorphous PE and 10.9 ± 1.9 xlO" 8 
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m 2 K/W for the oriented PE). They are also consistent with the results in Fig. 5 for the 
unfunctionalized tube for < 7 = 0. To examine this further, the radial distribution of the 
mass density, p, is shown in Fig. 6 for the three matrix types: the amorphous PE, the 
ordered PE and the EVA, each surrounding an unfunctionalized SWCNT. There is little 
distinction between the amorphous PE and the EVA. This is as expected, since there is 
relatively little vinyl acetate content in the EVA. The ordered PE simulation is clearly 
distinct from the other two, however. Thus, it is seen that modest changes in the matrix 
have little effect on the interfacial thermal resistance, at least for small ( 6 , 6 ) SWCNTs. 
Covalent functionalization of the tube, on the other hand, gave qualitatively different 
results, presumably by changing phonon coupling [3], 

3.3. Functionalized Graphite Nanoplatelet (GNP) Composites 

Fig. 7 shows calculated interfacial thermal resistances for both graphite surfaces 
as a function of the grafting density. One interesting feature is that the Rk values for the 
edge are substantially higher than for the flat graphite sheet surface, perhaps because the 
in-plane thermal motions cannot couple efficiently to the matrix. In contrast to results 
for a functionalized SWCNT [ 8 ] the length of the grafted chain appears to have no 
influence within the range studied here (n = 8-18). The net effect of increasing grafting 
density is also much smaller for the graphite surfaces than for the SWCNT. 

The interfacial thermal resistance values calculated with the MD simulations can 
be used with the effective medium equations of Section 2.3 to predict the effect of the 
functionalization on the thermal conductivity of the nanocomposite. The particles are 

o /: 

modeled as oblate ellipsoids with dimensions a 3 = 1x10" m and a\= 5x10" m. 
Considering first the functionalization of the larger flat surfaces of the GNP, Fig. 8 shows 
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thermal conductivity of the composite as a function of R„ , calculated using Eqs. 

" A perpendicular 

(2)-(8). Even with R K values approaching zero, the predicted thermal 

A perpendicular 

conductivity of the composite is only about 3 times that of the matrix at the highest 
loading of GNP (5 % volume). 

Considering next the effect of functionalizing the edges of the GNP, Fig. 9 shows 
predicted composite thermal conductivities as a function of R r . In this case, 

r r K parallel 

increases in thermal conductivity of over a factor of 15 are predicted at low values of 
R. r . The stronger dependence of thermal conductivity on /?„ versus 

R k ; ; is due to the greater conductivity of graphite along the parallel direction 

versus the perpendicular direction. The low R K a ; values in Fig 9. which would predict 

significant thermal conductivity increases are below those calculated by the MD 
simulation (indicated by the symbols) and may not be accessible experimentally. Thus, 
although functionalization of the faces and edges of GNPs can be expected to reduce 
interfacial resistances by 1/3 to 1/2, the effect on the thermal conductivity of the 
composite is not large. This result contrasts with the conductivity predictions for carbon 
nanotube composites in [8], which easily rose an order of magnitude above the base 
polymer. Both the thermal conductivity of the particle itself and the assumed aspect ratio 
were higher for the CNT. 

If the results for 0 and 20 weight % loading in reference [10] are interpolated to 
12 weight % ( = 5 volume %), the thermal conductivity obtained is 0.53 W/mK. This is 
in somewhat close agreement to the value plotted in Fig. 8 for 5 volume % and no 
functionalization (0.6 W/mK). Experimentally, the effectiveness of plate-like GNPs in 
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improving thermal conductivity varies considerably [10, 20, 21]. The analytical model is 
useful in exploring various parameters in addition to functionalization. For example, 
comparison of Figs. 9 and 10 shows the effect of changing the assumed aspect ratio of the 
ellipsoidal particle from 500 to 5000. All other values were unchanged. Although the 
composite thermal conductivities are substantially higher in Fig. 10, the overall 
improvement in thermal conductivity over the neat matrix is still modest. 

4. Conclusions 

Based on comparisons of amorphous and crystalline matrices, there does not 
appear to be any evidence that a “noncovalent functionalization” approach in the 
interfacial region will play a major role in reducing the interfacial thermal resistance and 
thereby increasing the thermal conductivity of the composite. 

MD results from EVOH functionalization of the SWCNT indicate that a variety of 
functionalization architectures are useful in achieving lower thermal interfacial resistance 
and therefore higher thermal conductivity. 

Simulations of GNP particles in EVA clearly predict that alkane functionalization 
should not result in an increase in thermal conductivity comparable to that seen with 
SWCNT filler. This is due to three factors: lower assumed conductivities of the GNP; the 
slightly higher interfacial resistance of the GNP (15-30 x 10" 8 m 2 K/W); and a smaller 
assumed aspect ratio. The multiscale approach presented here may be useful in 
optimizing nanocomposite properties. 
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Fig. 1. Views along the z-axis of atomistic models of the two polyethylene matrices 
surrounding a single wall carbon nanotube. The periodic boxes are 4.4 nm in the x and y 
dimensions and 5.7 nm in the (nanotube axial) z-axis direction. 
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Fig. 2. Modeling the flat graphite surface using periodic boundary conditions. The 
periodic boundary conditions are enacted across the x and z coordinates. The normal to 
the exposed surface is along the y axis. The periodic boundary dimensions are 4.7 nm in 
the x and z directions. Only the graphite atoms are shown in this depiction for clarity. 
The atoms associated with the functionalization and the matrix are not shown. 
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Fig. 3. Modeling the edge of the graphite particle using periodic boundary conditions. 
The periodic boundary conditions are enacted across the x and y coordinates. The normal 
to the exposed surface is along the z axis. The periodic boundary dimensions are 4.7 nm 
in the x and y directions. Only the graphite atoms are shown in this depiction for clarity. 
The atoms associated with the functionalization and the matrix are not shown. 
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Fig. 4 . The ellipsoid, depicted in two dimensions. The th ir d dimension, <22 = a\. 
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Fig. 5. The interfacial thermal resistance, Rk, as a function of the grafting density, a, of 
covalent chemical bonds attached to the wall of the carbon nanotube per unit area. 
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Fig. 6. The radial distribution of mass density, p, as a function of distance, r, from the 
center of the nanotube in atomistic nancomposite simulations. 
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Fig. 7. The interfacial thermal resistance as a function of the grafting density, <J, of 
covalent chemical bonds attached to the graphite nanoparticle. The open circles indicate 
values calculated for the edge surface of the parallel graphite sheets with a hydrocarbon 
chain of length n = 12 (GNP-C 12 H 25 ) covalently bonded. The filled diamond, square and 
triangle indicate values calculated for the flat surface of the graphite sheets with a 
hydrocarbon chain of length n = 8 , 12 and 18 (GNP-C ) ,H2„+i) covalently bonded, 
respectively. 
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Fig. 8 . The thermal conductivity of the GNP-polymer nanocomposite as a function of 
R„ for varying volume fraction,/. The dimensions of the GNP in this example 

A perpendicular 

O f\ 

are assumed to be a 3 = 1x10" m and a\ = 5x10" m. The dotted line indicates the value 
for the neat matrix (0.2 W/mK.) 
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R K parallel (xIOWK/W) 


Fig. 9. The thermal conductivity of the GNP-polymer nanocomposite as a function of 
R k for varying volume fraction,/. The dimensions of the GNP in this example are 

O /T 

assumed to be a 3 = 1x10" m and ci\ = 5x10" m. The dotted line indicates the value for 
the neat matrix (0.2 W/mK.) 
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R K parallel (xIOWK/W) 


Fig. 10. The thermal conductivity of the GNP-polymer nanocomposite as a function of 
R k n f for varying volume fraction,/. The dimensions of the GNP in this example are 

assumed to be a 3 = lxlO " 8 m and ci\ = 5xl0 " 5 m. The dotted line indicates the value for 
the neat matrix (0.2 W/mK.) 
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